PVA of the NBI: Study area

In this script we create the map for the study area.

Setup

## packages
library("tidyverse")
library("raster")
library("sf")
library("RStoolbox")
library("ggsn")
library("lwgeom")
library("rnaturalearth")
library("rworldmap")
library("rmapshaper")
library("shadowtext")
library("patchwork")
library("pdftools")
library("systemfonts")

## set the theme
theme_set(theme_minimal(base_family = "Open Sans"))
theme_update(
  axis.ticks = element_line(color = "grey10"),
  axis.ticks.length = unit(-.5, "lines"),
  axis.text.x = element_text(size = 16, 
                             margin = margin(t = 3)),
  axis.text.y = element_text(size = 16, 
                             margin = margin(r = 3)),
  panel.border = element_rect(color = "grey10", fill = NA, size = 1), 
  panel.background = element_rect(color = "white", fill = "white"),
  plot.margin = margin(5, 10, 10, 10),
  plot.tag = element_text(size = 24, face = "bold"),
  legend.position = "top",
  legend.text = element_text(face = "bold",
                             size = 18,
                             margin = margin(r = 15, b = 0))
)

Data

## elevation data
path_dem <- here::here("output", "geo-proc", "df_dem_alps.Rds")
path_hills <- here::here("output", "geo-proc", "df_hills.Rds")
path_slope <- here::here("output", "geo-proc", "df_slope.Rds")

if(!file.exists(path_hills)) {
  dem <- raster(here::here("data-raw", "geo-raw", "dem_3035_crop.tif"))
  dem_alps <- aggregate(dem, fact = 4)
  
  saveRDS(dem_alps, path_dem)
} else {
  dem_alps <- readRDS(path_dem)
}

## breeding and wintering locations
locations <- 
  tibble(
    location = c("B", "U", "K", "R", "L"),
    type = c("Breeding ground", "Breeding ground", "Breeding ground", "Breeding ground", "Wintering ground"),
    label = c("Burghausen\n(Germany)", "Überlingen\n(Germany)", 
              "Kuchl\n(Austria)", "Rosegg\n(Austria)", "Laguna di\nOrbetello\n(Italy)"),
    side = c("r", "l", "r", "r", "r"),
    x = c(4531605, 4258768, 4557635, 4629096, 4419572),
    y = c(2788488, 2739526, 2728565, 2616771, 2148807),
    from_x = c(4531605, 4258768, 4557635, 4629096, NA_real_),
    to_x = c(4419572, 4419572, 4419572, 4419572, NA_real_),
    from_y = c(2788488, 2739526, 2739526, 2616771, NA_real_),
    to_y = c(2148807, 2148807, 2148807, 2148807, NA_real_)
  ) 

Overview Map

sf_world <- 
  st_as_sf(rworldmap::getMap(resolution = "low")) %>% 
  st_transform(crs = st_crs(dem_alps)) %>% 
  st_buffer(dist = 0) %>% 
  dplyr::select(ISO_A2, SOVEREIGNT, LON, continent) %>% 
  mutate(area = st_area(.))

overview <- 
  ggplot(sf_world) +
    geom_sf(fill = "grey60",
            color = "grey80",
            lwd = 0.3) +
    geom_rect(xmin = 4000600, xmax = 4840000,
              ymin = 2000600, ymax = 2990000,
              color = "#212121",
              fill = NA) +
    geom_sf_text(data = sf_world %>% 
                   filter(ISO_A2 %in% c("DE", "IT", "AT")) %>% 
                   group_by(ISO_A2) %>% slice(1),
                 aes(label = ISO_A2),
                 family = "Open Sans",
                 color = "white",
                 fontface = "bold",
                 size = 5.7,
                 nudge_x = 20000) +
    geom_point(data = locations,
               aes(x, y),
               shape = 21,
               color = "#212121",
               fill = "white",
               size = 3,
               stroke = .8) +
    ggspatial::annotation_scale(location = 'tl',
                                text_family = "Open Sans",
                                text_cex = 1.2) +
    coord_sf(xlim = c(2000000, 6600000), 
             ylim = c(1000000, 5600000)) +
    scale_x_continuous(expand = c(0, 0), breaks = seq(-10, 30, by = 10)) +
    theme(panel.ontop = F,
          panel.grid.major = element_line(color = "grey85", linetype = "dashed", size = .4)) +
    labs(x = NULL, y = NULL)

overview

## save map
# ggsave(here::here("plots", "01_description", "NBIPVA_A_europe_countries.pdf"),
#        width = 8.95, height = 10, device = cairo_pdf)
# 
# pdf_convert(pdf = here::here("plots", "01_description", "NBIPVA_A_europe_countries.pdf"), 
#             format = "png", dpi = 750)

Hillshade map

dem_alps@data@values <- dem_alps@data@values * 10

if(!file.exists(path_hills)) {
  ## hillshade and slopeshade
  slope <- terrain(dem_alps, opt = "slope", unit = "radians")
  aspect <- terrain(dem_alps, opt = "aspect", unit = "radians")
  hill <- hillShade(slope, aspect, 40, 270)
  
  ## ggplot makes me turn the raster into points
  df_hills <- rasterToPoints(hill)
  df_hills <- data.frame(df_hills)
  colnames(df_hills) <- c("lon", "lat", "hills")
  
  ## the key to a pretty map is merging the slope shade with the hillshade
  df_slope <- rasterToPoints(slope)
  df_slope <-  data.frame(df_slope)
  colnames(df_slope) <- c("lon", "lat", "slope")
  df_slope$slope <- 1 - df_slope$slope ## invert the scale so that more slope is darker
  
  saveRDS(df_hills, path_hills)
  saveRDS(df_slope, path_slope)
} else {
  df_hills <- readRDS(path_hills)
  df_slope <- readRDS(path_slope)
}

## naturalearth data
oceans_raw <- 
  rnaturalearth::ne_download(scale = 10, 
                             category = "physical", 
                             type = "ocean",
                             returnclass = "sf")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\DataVizard\AppData\Local\Temp\Rtmp8AcCEd", layer: "ne_10m_ocean"
## with 1 features
## It has 3 fields
bbox <- as(extent(3, 17, 38, 51), 'SpatialPolygons')
crs(bbox) <- crs(oceans_raw)

oceans <-
  oceans_raw %>% 
  st_transform(crs = st_crs(bbox)) %>% 
  st_make_valid() %>% 
  st_crop(., bbox) %>% 
  ## transformation leads to inveresed oceans - fix it with ms_erase
  rmapshaper::ms_erase(st_as_sf(bbox), .) %>%
  st_transform(crs = st_crs(dem_alps))

lakes <- 
  rnaturalearth::ne_download(scale = 10, 
                             category = "physical", 
                             type = "lakes",
                             returnclass = "sf") %>% 
  st_transform(crs = st_crs(dem_alps)) 
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\DataVizard\AppData\Local\Temp\Rtmp8AcCEd", layer: "ne_10m_lakes"
## with 1355 features
## It has 41 fields
## Integer64 fields read as strings:  scalerank ne_id
## map
alps <- 
  ggplot(df_hills) +
    geom_raster(data = df_hills, aes(lon, lat, fill = hills, group = 1)) +
    geom_raster(data = df_slope, aes(lon, lat, fill = slope, group = 2), alpha = .6) +
    geom_sf(data = oceans,
            fill = "#627696",
            color = "#627696") +
    geom_sf(data = lakes,
            linetype = "dashed",
            fill = "#627696",
            color = "#627696") +
    # geom_curve(data = filter(locations, location %in% c("B", "K", "R")),
    #            aes(x = from_x, xend = to_x,
    #                y = from_y, yend = to_y),
    #            curvature = .12,
    #            size = 1.8,
    #            color = "#fefefe") +
    # geom_curve(data = filter(locations, location == "U"),
    #            aes(x = from_x, xend = to_x,
    #                y = from_y, yend = to_y),
    #            curvature = -.12,
    #            size = 1.8,
    #            color = "#fefefe") +
    geom_curve(data = filter(locations, location %in% c("B", "K", "R")),
               aes(x = from_x, xend = to_x,
                   y = from_y, yend = to_y),
               curvature = .12,
               size = .8,
               color = "#212121") +
    geom_curve(data = filter(locations, location == "U"),
               aes(x = from_x, xend = to_x,
                   y = from_y, yend = to_y),
               curvature = -.12,
               size = .8,
               color = "#212121") +
    ## outer ring
    geom_point(data = locations,
               aes(x, y,
                   color = type),
               shape = 21,
               fill = "#fefefe",
               size = 7,
               stroke = .8) +
    ## inner point
    geom_point(data = locations,
               aes(x, y,
                   color = type),
               size = 3) +
    ## location labels
    geom_shadowtext(data = filter(locations, side == "r"),
                    aes(x, y, label = label),
                    nudge_x = 20000,
                    size = 4.4,
                    color = "#fefefe",
                    family = "Open Sans",
                    fontface = "bold",
                    lineheight = .8,
                    hjust = 0,
                    bg.color = "#212121") +
    geom_shadowtext(data = filter(locations, side == "l"),
                    aes(x, y, label = label),
                    nudge_x = -20000,
                    size = 4.4,
                    color = "#fefefe",
                    family = "Open Sans",
                    fontface = "bold",
                    lineheight = .8,
                    hjust = 1,
                    bg.color = "#212121") +
    ggspatial::annotation_scale(location = 'tl',
                                text_family = "Open Sans",
                                text_cex = 1.2) +
    scale_color_manual(values = c("#d10173", "#21bffe"),
                       name = "") +
    rcartocolor::scale_fill_carto_c(palette = "ag_GrnYl", direction = -1, guide = "none") +
    #rcartocolor::scale_fill_carto_c(palette = "Fall", direction = -1, guide = "none") +
    scale_x_continuous(expand = c(0, 0), limits = c(4000600, 4840000)) +
    scale_y_continuous(expand = c(0, 0), limits = c(2000600, 2990000)) +
    guides(color = guide_legend(title.position = "left", 
                                title.hjust = 0.5, 
                                nrow = 1,
                                label.position = "right")) +
    theme(legend.position = "top") +
    labs(x = NULL, y = NULL)

alps

## save map
# ggsave(here::here("plots", "01_description", "NBIPVA_B_alps_GrnYl.pdf"),
#        width = 8.5, height = 10, device = cairo_pdf)
# 
# pdf_convert(pdf = here::here("plots", "01_description", "NBIPVA_B_alps_GrnYl.pdf"), 
#             format = "png", dpi = 750)
p <- overview + alps + plot_annotation(tag_level = "A", tag_suffix = ".")

p

## save panel
file <- here::here("plots", "01_description", "NBIPVA_panel_GrnYl")

ggsave(paste0(file, ".pdf"), width = 17, height = 10, device = cairo_pdf)

ggsave(paste0(file, ".svg"), width = 17, height = 10)

pdf_convert(pdf = paste0(file, ".pdf"), filename = paste0(file, ".png"), format = "png", dpi = 750)
## Converting page 1 to C:/Users/DataVizard/Google Drive/Work/Eco/Projects/2020_Drenske_PVA_NBI/plots/01_description/NBIPVA_panel_GrnYl.png... done!
## [1] "C:/Users/DataVizard/Google Drive/Work/Eco/Projects/2020_Drenske_PVA_NBI/plots/01_description/NBIPVA_panel_GrnYl.png"

Session Info
Sys.time()
## [1] "2022-06-30 17:57:47 CEST"
git2r::repository()
## Local:    master C:/Users/DataVizard/Google Drive/Work/Eco/Projects/2020_Drenske_PVA_NBI
## Remote:   master @ origin (https://github.com/EcoDynIZW/Drenske_2020_PVA_NBI.git)
## Head:     [af61606] 2022-06-28: update scripts
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
## [5] LC_TIME=German_Germany.1252    
## system code page: 65001
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] systemfonts_1.0.3   pdftools_3.0.1      patchwork_1.1.1    
##  [4] shadowtext_0.0.9    rmapshaper_0.4.5    rworldmap_1.3-6    
##  [7] rnaturalearth_0.1.0 lwgeom_0.2-8        ggsn_0.5.0         
## [10] RStoolbox_0.2.6     sf_1.0-5            raster_3.5-11      
## [13] sp_1.4-6            forcats_0.5.1       stringr_1.4.0      
## [16] dplyr_1.0.7         purrr_0.3.4         readr_2.0.2        
## [19] tidyr_1.1.4         tibble_3.1.6        ggplot2_3.3.6      
## [22] tidyverse_1.3.1    
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1         backports_1.2.1      spam_2.7-0          
##   [4] plyr_1.8.6           lazyeval_0.2.2       splines_4.1.0       
##   [7] jqr_1.2.1            listenv_0.8.0        geojsonlint_0.4.0   
##  [10] digest_0.6.29        foreach_1.5.1        htmltools_0.5.2     
##  [13] viridis_0.6.2        fansi_0.5.0          magrittr_2.0.1      
##  [16] geojsonsf_2.0.1      doParallel_1.0.16    tzdb_0.1.2          
##  [19] recipes_0.1.17       globals_0.14.0       modelr_0.1.8        
##  [22] gower_0.2.2          svglite_2.0.0        sysfonts_0.8.5      
##  [25] askpass_1.1          rmdformats_1.0.3     jpeg_0.1-9          
##  [28] colorspace_2.0-2     rvest_1.0.2          jsonvalidate_1.3.2  
##  [31] textshaping_0.3.6    haven_2.4.3          xfun_0.31           
##  [34] rgdal_1.5-28         crayon_1.5.1         jsonlite_1.7.2      
##  [37] survival_3.2-11      iterators_1.0.13     glue_1.4.2          
##  [40] gtable_0.3.0         ipred_0.9-12         V8_3.4.2            
##  [43] future.apply_1.8.1   maps_3.4.0           scales_1.2.0        
##  [46] qpdf_1.1             DBI_1.1.2            Rcpp_1.0.7          
##  [49] showtextdb_3.0       viridisLite_0.4.0    units_0.7-2         
##  [52] foreign_0.8-81       proxy_0.4-26         dotCall64_1.0-1     
##  [55] stats4_4.1.0         lava_1.6.10          prodlim_2019.11.13  
##  [58] httr_1.4.2           wk_0.6.0             geosphere_1.5-14    
##  [61] ellipsis_0.3.2       farver_2.1.0         pkgconfig_2.0.3     
##  [64] XML_3.99-0.8         nnet_7.3-16          sass_0.4.0          
##  [67] dbplyr_2.1.1         here_1.0.1           crul_1.2.0          
##  [70] utf8_1.2.2           caret_6.0-90         ggspatial_1.1.5     
##  [73] tidyselect_1.1.2     rlang_1.0.2          reshape2_1.4.4      
##  [76] munsell_0.5.0        cellranger_1.1.0     tools_4.1.0         
##  [79] cli_3.1.0            generics_0.1.2       ggmap_3.0.0         
##  [82] broom_0.8.0          evaluate_0.15        fastmap_1.1.0       
##  [85] ragg_1.1.3           yaml_2.2.1           ModelMetrics_1.2.2.2
##  [88] knitr_1.39           fs_1.5.0             s2_1.0.7            
##  [91] RgoogleMaps_1.4.5.3  showtext_0.9-4       future_1.22.1       
##  [94] nlme_3.1-152         xml2_1.3.2           compiler_4.1.0      
##  [97] rstudioapi_0.13      curl_4.3.2           png_0.1-7           
## [100] e1071_1.7-9          reprex_2.0.1         bslib_0.3.1         
## [103] stringi_1.7.5        highr_0.9            fields_13.3         
## [106] rgeos_0.5-9          lattice_0.20-44      Matrix_1.3-3        
## [109] classInt_0.4-3       vctrs_0.3.8          pillar_1.7.0        
## [112] lifecycle_1.0.1      jquerylib_0.1.4      geojsonio_0.9.4     
## [115] data.table_1.14.2    bitops_1.0-7         maptools_1.1-2      
## [118] rcartocolor_2.0.0    R6_2.5.1             bookdown_0.24       
## [121] KernSmooth_2.23-20   gridExtra_2.3        parallelly_1.29.0   
## [124] codetools_0.2-18     MASS_7.3-54          assertthat_0.2.1    
## [127] rprojroot_2.0.2      rjson_0.2.21         httpcode_0.3.0      
## [130] withr_2.5.0          geojson_0.3.4        parallel_4.1.0      
## [133] hms_1.1.1            terra_1.4-22         rpart_4.1-15        
## [136] timeDate_3043.102    class_7.3-19         rmarkdown_2.11      
## [139] git2r_0.29.0         pROC_1.18.0          lubridate_1.8.0